module double_kind
  implicit none
  integer, parameter :: dp = 8
  integer, parameter :: nODE = 2  
end module double_kind

module inteface_ydot
  interface
     function ydot(t, w) result(y)
         use double_kind
         implicit none
         real(dp), intent(in) :: t
         real(dp), intent(in), dimension(nODE) :: w
         real(dp), dimension(nODE) :: y
       end function ydot
    end interface
  end module inteface_ydot

module interface_rk4
  interface
    function rk4(dt, w0, t) result(rsl)
       use double_kind
       use inteface_
       implicit none
       real(dp), intent(in) :: dt, t
       real(dp), intent(in), dimension(nODE) :: w0
       real(dp), dimension(nODE) :: rsl, s1, s2, s3, s4
     end function rk4

   subroutine print_matrix(A)
      integer, parameter :: dbl = 8
      real(dbl), dimension(:,:), intent(in) :: A
    end subroutine print_matrix
    
  end interface
end module interface_rk4
  
       
  

     
program rk4_2
  use double_kind
  use interface_rk4
  implicit none
  integer, parameter :: step = 10
  real(dp), dimension(step+1, nODE) :: y
  real(dp), dimension(step+1) :: t

  real(dp) :: dt
  integer :: k

  t(1) = 0
  dt = 1.0/step
  y(1,1) = 0.0
  y(1,2) = 1.0
  do k = 1, step
     t(k+1) = t(k) + dt
     y(k+1,:) = rk4(dt, y(k,:), t(k))
  end do

  call print_matrix(y)



  
end program rk4_2

function rk4(dt, w0, t) result(rsl)
  use double_kind
  use inteface_
  implicit none
  real(dp), intent(in) :: dt, t
  real(dp), intent(in), dimension(nODE) :: w0
  real(dp), dimension(nODE) :: rsl, s1, s2, s3, s4
 ! real(dp), external :: ydot  
  
  s1 = ydot(t, w0)
  s2 = ydot(t+dt/2, w0+dt*s1/2)
  s3 = ydot(t+dt/2, w0+dt*s2/2)
  s4 = ydot(t+dt,   w0+dt*s3  )

  rsl = w0 + dt*(s1 + 2*s2 + 2*s3 + s4)/6
end function rk4

function ydot(t, w) result(y)
  use double_kind
  implicit none
  real(dp), intent(in) :: t
  real(dp), intent(in), dimension(nODE) :: w
  real(dp), dimension(nODE) :: y

  y(1) = w(2)*w(2) - 2*w(1)
  y(2) = w(1) - w(2) - t*w(2)*w(2)
end function ydot

subroutine print_matrix(A)
  integer, parameter :: dbl = 8
  real(dbl), dimension(:,:), intent(in) :: A
  integer :: row, k
  row = size(A,1)
  do k = 1, row
     write(*,*) A(k,:)
  end do
  write(*,*)
end subroutine print_matrix
